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Abstract 

The new version of the HZETRN deterministic transport code based on Green’s function 
methods, and the incorporation of ground-based laboratory boundary conditions, has lead 
to the development of analytical and numerical procedures to include off-axis dispersion 
of primary ion beams due to small-angle multiple Coulomb scattering. In this paper we 
present the theoretical formulation and computational procedures to compute ion beam 
broadening and a methodology towards achieving a self-consistent approach to coupling 
multiple scattering interactions with ionization energy loss and straggling. Our initial 
benchmark case is a 60 MeV proton beam on muscle tissue, for which we can compare 
various attributes of beam broadening with Monte Carlo simulations reported in the open 
literature. 

1. Introduction 

Current developments in the High Charge and Energy Transport (HZETRN) model are 
focused towards a full three-dimensional and computationally efficient deterministic 
transport code, capable of simulating HZE radiation transport with either space or 
ground-based laboratory boundary conditions (Wilson et ah, 2005). Two parallel 
development paths of HZETRN have ensued. The first development path for HZETRN 
was based on marching procedures adequate for boundary conditions found in the space 
environment. For nearly two decades, this version of HZETRN has been validated in the 
space environment by comparing to a wide variety of spaceflight measurements taken on 
the Space Transportation System (STS) and the International Space Station (ISS). A 
complementary development path for HZETRN has been initiated in recent years based 
on a Green’s function approach (Tweed et ah, 2005) combined with nonperturbative 
techniques (Wilson et ah, 1994). In addition to computation efficiency, an advantage of 
the Green’s function approach is the ease and flexibility by which both ground-based 
laboratory and space environment boundary conditions can be incorporated into a single 
theoretical formulation. Ground-based laboratory boundary conditions are characterized 
by directed beams of ions of specific energy with detailed diagnostic particle 
spectrometer devices. 

One aspect of the new version of HZETRN is the incorporation of off-axis dispersion of 
the primary ion beam for ground-based directed-beam applications. The beam broadening 
mechanism is attributed to small-angle, multiple Coulomb scattering of the incident ions 
by the target material nuclei. While the effects of multiple scattering are negligible in the 
space radiation environment, multiple scattering must be included in directed-beam 
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applications and ground-based laboratory transport code simulations. To name a few 
applications, multiple scattering mechanisms must be included in directed-beam ion 
transport codes in order to accurately model and interpret the response of particle 
spectrometer devices in ion beam experiments, to simulate the physical and biologically 
important radiation dose, and to develop new methods and strategies for light ion beam 
radiation therapy. 

In this paper we present the theoretical formalism and computation procedures for 
incorporating small-angle, multiple Coulomb scattering into the deterministic HZETRN 
code, and the method of coupling the ion-nuclear scattering interactions with ionization 
energy loss and straggling. Simulations of the effects of multiple scattering on ion beam 
broadening characterization will be computed for a 60 MeV proton beam incident on 
muscle tissue. This example has important practical applications for improving 
knowledge in radiation biology and proton beam radiation therapy. Furthermore, this 
example will provide our initial benchmark to compare our new computation procedures 
with Monte Carlo simulations reported by Noshad and Givechi (2005). We calculate the 
effects of multiple scattering on path-length corrections, lateral beam broadening, and 
absorbed dose. 

2. The Boltzmann Transport Equation 

Beam broadening is attributed to the accumulation of many random, small-angle 
scattering events by which a particle from the incident ion beam is deflected by the 
average Coulomb potential established by the constituent nuclei of the target material and 
partially shielded by the orbital electrons. Thus, the propagation and broadening of an ion 
beam due to a random distribution of many scattering events, as described above, can be 
formulated in terms of the transport of a statistical distribution function. Specifically, 
transport of HZE ions is described by the linear Boltzmann equation, which can be 
derived on the basis of conservation principles (Wilson et ah, 1991). Considering only 
ion-nucleus Coulomb scattering interactions, the transport of an ion beam is described by 

U • V3>(x,tJ,£) =jCT(U,U',E,E / )4)(x,U',£'Vu'^'-cr(£')<t>(x,U,£’), (1) 

where <E(x,U,£') is the flux of ions at position x, propagating in direction Q with energy 

E. The differential and total macroscopic cross sections for ion-nucleus Coulomb 
scattering are denoted cr(U,U',E,E') and o(E) , respectively. The first term on the right- 

hand-side of (f) describes the production per volume of particles at position 
x propagating in direction U with energy E due to Coulomb scattering of particles 
initially at position x propagating in direction U' with energy E ' . The second term on the 
right-hand-side of (1) describes of the loss of particles per volume due to scattering of 
particles at position x in direction U with E out of the beam. The net gain of particles per 
volume due to scattering, as described by the right-hand- side of (1), is identified with the 
infinitesimal change per unit length of particle flux at position x propagating in direction 
U with energy E, which is described by the left-hand-side of (1). 

3. Reduction to a Fokker-Planck Diffusion Equation 
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In this section we invoke two approximations that simplify both the scattering cross 
section and the Boltzmann transport equation in (1). The first approximation is that the 
Coulomb scattering is elastic. Thus, the differential Coulomb scattering in (1) takes the 
form 


o(t,V',E,E') = o(t,t',E)8(E-E'). 


( 2 ) 


The differential scattering cross section for elastic Coulomb scattering from a fixed target 
nuclei is the Rutherford differential scattering cross section (Rossi and Greisen, 1941). In 
the small-angle scattering limit, our second approximation, the Rutherford differential 
cross section depends only on the scattering angle - i.e., the angle between the incident 
particle direction and the direction of the scattered particle, and is given by 


a(U,U',C) = a(U«U',^) = cr(0,C) = 4^ 


a t (v- P y o 4 


( 3 ) 


where Na is Avogadro’s number; A T is the target molecular weight; Z T and Z P are the 
target and projectile atomic numbers, respectively; r 0 is the classical electron radius; u e is 
the electron rest mass; v and p are the projectile velocity and momentum, respectively; 
and 0 is the scattering angle. Note that the expression for the differential cross section 
implies that our unit of length is g/cm 2 . 

The elastic scattering approximation and the small-angle limit reduce the form of the 
differential scattering cross section to the expressions in (2) and (3), which simply the 
transport equation in (1). The small-angle limit will yield additional simplifications to the 
transport equation, as outlined below. 

Consider the ion beam incident at the origin and propagating along the positive z-axis. 
The components of the unit vector oriented along the propagating direction, in terms of 
our coordinate system, are given by 

U = sin0cos/3v + sin0sin/3v + cos0z (4) 

where 0 is the angle between the z-axis and Q, and (3 is the azimuth angle. It’s useful to 
define orthogonal projections of the scattering angle (0) on a plane perpendicular to the 
initial beam direction (i.e., the z-axis): 

d x = sin 9 cos /3 
0 = sin 0 sin (>. 

In the small-angle limit of d 1 , the projection angles are approximated by 

0 v =0cos/3 (7) 


( 5 ) 

( 6 ) 
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8 y « d sin /3 


( 8 ) 


such that 


e 2 ~e 2 x+ d 2 y . 


(9) 


Moreover, the drift operator on the left-hand side of (1) reduces to 


U 


•v~e x —+e y —+—. 

dx ' dy dz 


( 10 ) 


Considering the approximations and definitions discussed above - i.e, the fomi of the 
differential elastic scattering cross section in (2) and (3) and the drift operator in (10) - 
the transport equation in (1) can be written as 


a a a 
- + e x - + e v — 

dz dx ' dy 


®(x,y,z;6 x ,d y ) d/3 £ o(d)ddd ^®(x,y,z;d x ,d y ')-®(x,y,z;d x ,d y ) j (11) 


A diffusion equation can be derived by making the usual Fokker-Planck approximation 
by expanding the expression in brackets, (1 1), in a Taylor series and retaining no terms 
beyond the second derivatives (Scott, 1949). The Fokker-Planck approximation is 

justified in the small-angle limit. Realizing that (8 X ,6 X ) represent the projected angles 

before scattering into ( 6 x , 6 y ) we have 

0,'= 0,-0 cos p (12) 

0,'- 0,-0 sin 13. (13) 


Thus, the Taylor series expansion of the ion flux distribution is written 

t t <30 <30 

0(x, y, z; 8 x , d y ) = 0(x, y, z; 0 v , Q y ) - 0 cos p — - 0 sm p — 

dU x 66 v 

0 2 cos 2 /3d 2 O d 2 sin 2 6 0 2 O 
-| ! 1 1 

2 3d 2 2 68 2 

x y 

Substitute (14) into the right hand side of (1 1) and integrate over solid angle 
(dQ' » dddd/3 ). The result is 


(14) 


d d d 

, 1 

' d 2 d 2 

h 8 X — + 9 y — 

dz dx ' dy 

O (x,y,z;d x ,6 v ) = 

A s (z) 

_dd 2+ dd 2 _ 


O (x,y,z;6 x ,d y ) 


(15) 
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where 




1 rf(e 2 U)) 

4 dz 


= 16jt 7V, (Z/Z/,) (/ ” jU< ) , ln(181Z: 1/3 ). 

4 (v(z)^(z)) 2 


( 16 ) 


The characteristic scattering length is given by k s , which is equal to the inverse of the 
diffusion coefficient, and is related to the mean- square scattering angle per unit thickness 
(dz), as (16) indicates. Equation (16) can be readily extended to include multi-component 
target materials, such as muscle tissue discussed below, by the expression 

= \6nZ 2 p , y P A (m)Z 2 (m) ln(1 8 1Z~ 13 (m)) (17) 

K(z) (v(z)^(z)) V 


where the sum is over the chemical constituents of the target material and Pa is the 
number of atoms per gram. 

We have tacitly coupled energy loss in the Fokker-Planck partial differential equation in 
(15) by letting the diffusion coefficient depend on the propagation depth (z), and by 
letting the scattering angle depend on propagation depth in (16). The expression for the 
characteristic scattering length in (16) was evaluated using the approximations employed 
by Rossi and Greisen (1941): the lower limit of integration in (16) was set by considering 
the screening of the nuclear charge by the orbital electrons using the Thomas-Fermi 
statistical model of the average Coulomb potential function; and the upper limit of 
integration in (16) was set by considering the finite nuclear size and assuming a uniform 
spherically symmetric nuclear charge distribution. 

The energy loss of a projectile particle is due to energy transfer to the orbital electrons of 
the target material. Using the continuous slowing down approximation, the mean energy 
of the ion beam after propagating a distance z into the target material is given by the 
usual range-energy relation (Wilson et al., 1991) 

{E{z)) = R-'[R{E 0 )-z\ (18) 


where R(Eo) is the mean range of an incident particle at the initial energy Eo at zero 
depth, and the unit of energy is taken to be MeV/nucleon. The velocity and momentum at 
each propagation depth-z in (16) can be evaluated using the mean beam energy from (18) 
and the relativistic energy relations, such that 


/3(z) = v(z)/c = 


1 - (l + (C (z)) / amu • c 2 ) 


(19) 
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p(z)c = A p • (amu • c 2 ) ■ p (z) • (l - /3 2 (z) ) ' “ . 


( 20 ) 


In the above equations, c is the speed of light, A p is the molecular weight of the projectile 
particle, and amu = 931.5 Mev/c 2 (atomic mass unit). 

4. Solution of Depth-Dependent Fokker-Planck Equation 

The 3-D Fokker-Planck equation in (15) can be separated into two independent partial 
differential equations in the two sets of projected coordinates - x, 0x and y, 0y - by 
factoring the flux distribution into a product of projected flux distributions, i.e., 

®(x,y,z;6 x ,6 y ) = O(z,x;0 v )O(z,y;0 v ). (21) 

For the simple boundary conditions of an infinitesimal pencil beam of unit flux at the 
origin, i.e., 


d>(O,x,0 v ) = <5(x)<5(0 A ) (22) 

<D(0 ,y,0 y ) = d(y)6(0 y ), (23) 


Eyges ( 1 948) found the solutions of the depth-dependent Fokker-Planck partial 
differential equations, which are referred to as the Fermi-Eyges distribution functions. 
Thus, the solutions for the projected flux distributions are 


<h(z,x,0 T ) = 


<E(z,y,0 v ) = 




y’ 4 njB{z) 


exp 


exp 


~(A 0 (z)x 2 -2xd x A ] (z)+A 2 (z)0}) 
4 B(z) 

-(A o (z)y 2 -2y0 y A 1 (z)+A 2 (z)6 2 ) 
4 B(z) 


(24) 

(25) 


where 


B(z) = A 0 (z)A 2 (z) - (z). 


(26) 


The projected- flux distributions in (24) and (25) resemble 2-variate Gaussian probability 
distributions in angular and lateral displacements, with a non-zero correlation between 
the angular and lateral displacements. The similarity of the projected- flux distributions 
with 2-variate Gaussian probability distributions is consistent with our physical concept 
and formulation of a localized initial ion beam that is broadened with propagation 
distance in its angular width and lateral displacements due to a continuous distribution of 
random, small-angle scattering events. This analogy is further strengthened by the fact 
that the A-functions in (24) and (25) are integral moments of the diffusion coefficient, 
which can be expressed in terms of the usual statistical moments. Specifically, 
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(27) 


24M -(<?«> - 2 filfc-if' 

2A l (z).{ye,(z)).2j i a 


z dld 2 (z') 


dz' 


dz' 


(z-z')dz' _ i 
10 As( z 9 


d(e 2 (z')\ 

= jJ q (z ~ z) dz 


2 4 (z).(/( Z) >=2j;<^-ij;Vz') 


dz’ 


d/e 2 (z')) 

'^ 2 [ 1 dz'. 


dz' 


(28) 

(29) 


Similar equations hold for the statistical moments projected on the z-x and z-0x planes. 

To accommodate realistic boundary conditions, and in order to calculate the path-length 
corrections in section 5, more general solutions for the Fermi-Eyges distributions are 
required. We extend the boundary condition in (22) and (23) to describe an infinitesimal 
pencil beam of unit flux with both non-zero angular displacements (0x o and 0y o ) and non- 
zero lateral displacements (x 0 andy 0 ). That is, we take 


d>(0, x, 6 X ) = d(x- x 0 )8 (6 X - 6 Xo ) (30) 

^(O, y, d y ) = d(y- v 0 )<5 (d y - 0 Vq ). (31) 

Generalized Fermi-Eyges projected- flux distribution functions can be obtained by solving 
the projected Fokker-Planck partial differential equations subject to the boundary 
conditions expressed in (30) and (31). The form of the generalized Fermi-Eyges 
distribution functions are the same as (24) and (25) provided the following variable 
replacements are made: 


*-*•(*- *o - e x 0 z ); -*• ( & x ~ d x 0 ) (32) 

y (y - To - e yo z); e y - (p y - e yo ). (33) 

The statistical moments in (27)-(29) can also be computed from the generalize Fermi- 
Eyges distribution functions. Statistical boundary conditions can be accommodated by 
replacing the initial angular and lateral displacements by an ensemble average such that 
(Hollmark et ah, 2004) 


(0 2 (z)} = (0 2 (O))+J^^!iiVj r ' 


dz’ 


(r 2 (z)) = (r 2 (0)) + 2 (rd( 0)) z + (d 2 (0)} z 2 +£ (z - z') 2 


dle 2 (z') 


(rd(z)) = (rd( 0)> + (0 2 (O))z + (z - z')- 


dld 2 { z') 


dz' 


dz’ 


dz' . 


dz' 


(34) 

(35) 

(36) 


The above equations, (34)-(36), are transport integrals that describe beam broadening 
characteristics for arbitrary boundary conditions. We have assumed, for simplicity, 
cylindrical symmetry in deriving (34)-(36). Asymmetries can be accounted for by 
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developing the transport integrals for the statistical moments in the projected z -x/0x and 
z-y/0y planes separately. For this study, however, we will assume cylindrical symmetry. 

In Figure 1 and Figure 2 we show, respectively, the radial broadening and the radial flux 
distribution for an initial infinitesimal 60 MeV proton beam incident on muscle tissue. 
The longitudinal propagation extends to 31 mm. The Bragg peak is located at 30.49 mm 
(shown in Figure 8). The chemical composition of muscle tissue was considered as 
H(63%)+C(6%)+N(l%)+0(28%), and its density was considered to be equal to 1 g/cm 3 
(Noshad and Givechi, 2005). The radial flux distribution is obtained from the projected 
Fermi-Eyges distributions functions by 


+oo +oo 

f(x,y,z ) =J dd xf^ d0 } M z ,x,d x )®(z,y,d y ) 
1 (~(x 2 +y 2 ) 


n(r 2 ( z )) P | (r 2 (z)) 


(37) 


Results for the 60 MeV proton beam on muscle tissue described above were computed by 
the TRIM Monte Carlo code and reported by Noshad and Givechi (2005). The radial 
broadening at the Bragg peak in Figure 1 is 2. 13 nun, which is in excellent agreement 
with the results calculated by Noshad and Givechi (hereafter referred to as NG). 

However, the broadening at intermediate tissue depths and the rate of broadening is 
substantially different between NG and our calculations. Our radial broadening is greater 
than the NG results for tissue depths less than the Bragg Peak. These differences are 
relatively small below 20 mm tissue depth, but appear to be on the order of 50% at 
around 25 mm. The radial broadening computed by NG rapidly increases between 25 mm 
and the Bragg Peak. 

Figure 2 shows the radial flux distribution for an initial infinitesimal 60 MeV proton 
beam on muscle tissue at various tissue depths. The flux is shown along one of the lateral 
axis (either x- or y-axis). The proton beam is highly localized at 5 mm depth where the 
radial broadening is only 0.12 mm. At 15 mm, 25 mm, and 30 mm tissue depths, the 
radial broadening is 0.64 mm, 1.48 mm, and 2.06 mm, respectively. 

In Figure 3 and Figure 4 we show, respectively, the radial broadening and the radial flux 
distribution for an initial finite width 60 MeV proton beam on muscle tissue. The initial 
beam width is 2 mm in diameter. In this case our radial broadening results are less than 
those computed by NG. The largest differences occur for tissue depths greater than 20 
mm. We find a radial broadening of 2.35 mm at the Bragg Peak, while the radial 
broadening is closer to 3.00 mm for the result reported by NG. Similar to the initial 
infinitesimal beam case, the NG radial broadening increases at a slow rate at small to 
intermediate tissue depths, but increases rapidly near the end-of-range, particularly 
between 20 mm and the Bragg peak for the initial 2 mm diameter proton beam. 

Figure 4 shows the radial flux distribution for an initial 2 mm diameter 60 MeV proton 
beam on muscle tissue at various tissue depths. The flux is shown along one of the lateral 
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axis (either x- or y-axis). The radial broadening at 5 mm, 15 mm, 25 mm, and 30 mm 
tissue depths is 1.01 mm, 1.19 mm, 1.79 mm, and 2.29 mm, respectively. 

5. Path-Length Correction 

To summarize the theoretical formulation of sections 3 and 4, energy loss was coupled to 
beam broadening by solving the multiple scattering transport equation in (15) with the 
depth-dependent diffusion coefficient given by (17), using (19)-(20). The diffusion 
coefficient is a function of depth due to ionization energy loss and is calculated assuming 
the continuous slowing down approximation in (18). However, the actual particle path in 
the target material is greater than the physical depth-z in the material due to multiple 
scattering processes. A necessary step towards a self-consistent theory of multiple 
scattering and energy loss is to compute a path-length correction to the longitudinal depth 
in the evaluation of the mean energy of the beam particles, prior to computing the 
transport integrals in (34)-(36) and the flux distributions in (24)-(25), and (37). 

Bichsel and Uehling (1960) defined an effective 1-D path-length correction which is a 
function of the projected mean-square scattering angle, 

(A R(Z)) =j^(d;(z'))dz’ =f o (d 2 y (z'))dz'. (38) 

The projected mean-square scattering angle at intermediate depths- z depends on the 
boundary conditions at zero depth and at depth-z, which can be formulated in terms of a 
doubly-conditional probability distribution 


(e 2 M)) =C p(d >' 1 °>y<» e o’>wA K d0 'y ( 39 ) 

The probability function in (39) describes the probability that at intermediate depth-z' a 
beam particle is propagating with an angular divergence of 6 ' v , given that the particle had 

lateral displacement y 0 and angular displacement 0 O at zero depth (z=0), and emerges at 
depth-z with a lateral displacement yi and an angular displacement 0i. 

The influence of initial beam statistics has already been included in the development of 
the transport integrals in (34)-(36). Therefore, without loss of generality, we can set y 0 
and 0 O equal to zero in (39). The emerging boundary conditions at depth-z require a little 
more insight. Recall that the goal is to develop an effective 1 -D propagation depth in the 
target material. Thus, we develop a correction to the physical depth-z - i.e., the effective 
1-D propagation depth - by computing a path-length correction due to multiple scattering 
of particles that emerge at depth-z with zero angular displacement, irrespective of the 
lateral displacement. Consequently, for the boundary conditions just described, the 
doubly-conditional probability distribution in (39) becomes 
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( 40 ) 


+oo 

p(e' y (z'))de' y d yi p(e' y -z'\o,o,o-z,y x ,o) 

+00 +00 

d6 'yJ^ d y\ f , dy'<l>(z - z',y l - y - d'.(z - z'),-d' y )<&{z ,y ,d ' y ) 

+oo 

f dy^(z,y v 0) 


= dQ' 

y 


( Mz) lv P 

(< 

[ 1 1 1 \ 
| 

\\x\{z-z)\(yz) j 

4 

\ 

Mz-z') M z ')\j 


Substitute (40) into (39), and (39) in to (38), such that 

2 A( Z ')A( Z ~ Z ') 


(«; 2 < z ')>- 


4 >( z ) 


(41) 


and 


(A R(z)) = 



Ag(z - z')A 0 (z')dz'. 


(42) 


To get a rough idea of how the path-length correction in (42) depends on depth and the 
characteristic scattering length in (17), assume the scattering length is independent of 
depth. The statistical moments in (27)-(29) can be evaluated analytically, and (42) 
reduces to 


(A R(z)) 



(43) 


Figure 5 shows the percent increase in the penetration depth due to small-angle multiple 
scattering. At the end-of-range, the effective 1-D path-length is 0.8% greater than the 
physical depth-z. This may appear as a small correction, but it has a significant effect on 
the mean energy, which is decreasing non-linearly near the end-of range. The influence 
of multiple scattering on the mean energy is characterized by replacing (18) with the 
expression 


(■ E(z)) = R - 1 [/?(£„) -(z + (AK(z)))]. 


(44) 


Figure 6 shows the influence of the path-length correction on the mean energy. Since the 
effective 1-D propagation depth is greater than the physical depth-z, the mean energy due 
to multiple scattering processes is less than the mean energy without scattering. 
Furthermore, the ion beam will reach the end-of-range at a smaller physical depth when 
the effects of multiple scattering on the propagation depths are included. For example, the 
minimum mean energy at the location of the Bragg peak when scattering effects are 
included is 0.14 MeV/nucleon, and the Bragg peak is located at 30.49 mm tissue depth. If 
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scattering effects are ignored, the mean energy at 30.49 mm is 3.95 MeV/nucleon. This 
corresponds to nearly a factor of 30 overestimation of the mean energy due to the meager 
0.8% increase in effective 1-D path-length due to multiple scattering processes, as 
described above. However, the overestimation in the Bragg peak location, if multiple 
scattering effects on path-length are ignored, is comparable to the error in effective 1-D 
path-length. Thus, the Bragg peak is located at 30.73 mm tissue depth if multiple 
scattering is ignored, corresponding to an overestimation of roughly 0.8% 

Energy straggling is also influenced by multiple scattering, through the effective 1-D 
path-length correction, since straggling is a function of the mean energy (Wilson et al. 
2000, 2002). Figure 7 shows the increases in energy straggling as a result of multiple 
scattering path-length corrections for the 60 MeV proton beam on muscle tissue. The 
energy straggling width increases with decreasing mean energy (Wilson et al., 2000, 
2002). A smaller mean energy when multiple scattering effects are includes produces a 
larger straggling width. At the Bragg peak when multiple scattering processes are 
included (i.e., at 30.49 mm tissue depth), the straggling width is 2.21 MeV. If multiple 
scattering effects are not included, the straggling width is 2.10 MeV at 30.49 mm tissue 
depth. 

6. Absorbed Dose Distribution 

In this section we describe the formulation to compute absorbed dose, which includes the 
coupling of multiple scattering with ionization energy loss and straggling. 

6. 1 On-Axis Dose Distribution 

The on-axis dose distribution is given by the following expression 


oo 

d(z) = Kf o S(E)(D(z,E)dE, 


(45) 


where S(E) is the stopping power, d>(z,E) is the on-axis spectral flux distribution, and K is 
a unit conversion factor (Wilson et al., 2005). The spectral flux distribution is obtained by 
a second-order moment expansion of the Boltzmann transport equation, which can be 
expressed in terms of a Gaussian function 


d>(z,£) 



2 s 2 (z) 


(46) 


where the mean energy in (46) is evaluated using (44) and the energy straggling width 
s(z) is evaluated according to the global Payne formulation described by Wilson et al. 
(2000, 2002). Both the mean energy and straggling width include multiple scattering 
effects via the path-length correction, as describe in the previous section. An expression 
for the absorbed dose can be obtained by substituting (46) into (45): 
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d(z) = 


K 




2\[n 


HS 


{■ E(z)) + y[2s(z)x h m ] 


(47) 


for penetration depths less than the particle range, and 
K 


d(R + A z) = 


2 yfn | 
/f 

2yfjl 


i E{z) ) V G' exp 
yf2s(z)^" n -° 



S 


8 s 2 (z) K m) 


2 V mj 


(48) 


X' ^ £7 O 

2jm=0^ m 


'{E(z)) + j2s(z)x h m } 


for penetration depths greater than the end-of-range. The above expressions assume 
incident unit particle flux. In (48) R is the actual mean range of the beam, including the 
effects of multiple scattering, and Az is the increment in depth beyond the mean range. In 
both (47) and (48), H m and x h m are the m th Gauss-Hermite quadrature weight and abscissa, 
respectively. In (48), G m and xf n are the m th Gauss-Legendre quadrature weight and 

abscissa, respectively. The terms in (48) that involve the Gauss-Legendre weights and 
abscissas account for the statistical distribution of particles with energies greater than 
zero when the mean energy reaches zero at the end-of-range. The number of Gauss- 
Hermite and Gauss-Legendre quadrature points are denoted qh and q g , respectively. 

Note: one can show that the orthogonal polynomials over the range [0,°°] are the standard 
Gauss-Hermite polynomials with the numerical quadrature formula specified by 
evaluating the standard Gauss-Hermite abscissas at only the positive values and to divide 
the standard quadrature weights by two. 

Figure 8 shows the relative absorbed dose for the 60 MeV proton beam on muscle tissue. 
The dose-depth curve in Figure 8 was calculated from (47) and (48) using a 5-point 
Gauss-Hermite quadrature rule and a 10-point Gauss-Legendre quadrature rule. The 
result is in excellent agreement with the Monte Carlo TRIM code calculation reported by 
NG. In section 5 it was shown by the minimum mean energy that the location of the 
Bragg peak is at a tissue depth of 30.49 mm, which is consistent with the maximum in 
absorbed dose in Figure 8. NG report the Bragg peak at 30.5 mm. Recall from section 5 
that the Bragg peak would be located at 30.73 mm tissue depth if multiple scattering 
effects via path-length corrections was not included in the calculations. 

6.2 Off-Axis Dose Distribution 

The conceptual picture of ion beam transport is that the beam is composed of a collection 
of independent pencil beams propagating and broadening through the target material. 
Thus, the dose distribution is a summation of contributions from each pencil beam. This 
concept of the off-axis dose distribution can be formulated as follows (Hogstrom et al, 
1981): 


mx,y,z).flS(x',y,z)d(x-x,y-y',z)dx’dy. 


(49) 
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In the above equation, S(x’, y') is the relative strength of a pencil beam at x and y . The 
contribution of dose at x and y from the pencil beam at x and y is specified by the dose 
function d(x - x',y - y',z ) . The integration limits in (49) are specified by the 2-D lateral 
extent of the beam projected at depth-z, which includes the initial beam width and 
broadening due to multiple scattering. 

The dose function in (49) is factored into an on-axis term, given by (45), which is 
modulated by an off-axis beam broadening function, given in terms of the radial Fermi- 
Eyges distribution function, such that 

d (x, y, z) = /(x, y, z)d (z) (50) 

where 


+oo +oo 

f(x,y,z) =J dd x J_^ dd y ®(z,x,d x )®(z,y,d y ) 
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no; (z) 


exp 


-(* 2 + y 2 ) 

<T 2 (z) 


(51) 


and 


o;(z) =J o (z-z) 




(52) 


Notice that the multiple scattering distribution function in (51) and the transport integral 
in (52) have slightly different forms compared to the previous expressions in (37) and 
(35). In computing the off-axis dose distribution as formulated in (49), the influence of 
the finite initial boundary conditions are included in the limits of integration rather than 
in the transport integrals. Thus, the beam broadening function and the transport integral 
in (51) and (52) need only include broadening due to multiple scattering as a result of 
transport through the target material. The boundary condition on the double integral in 
(49) is the projected lateral extent of the broadened beam at depth-z, as previously 
mentioned. 


As an example of an application of the above formulation, consider an initial finite beam 
of uniform unit flux distribution incident on the target material at z = 0. The relative 
strength of a pencil beam is specified by the following equation 

S(x,y,z) = 9 (x + 4 + )? (-* + f + V a r( Z ) y ( V + f + j°r( Z ) y(^y+i+ \l°r( Z ) ) ( 53 ) 


In the above equation, d is the diameter of the initial beam, and 0 in this case denotes the 
Heavyside step function. Substitute (53), (50), and (51) into (49) to obtain the following 
result 
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D(x,y,z) = ^ 
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where erf(x) is the error function, and the transport integrals are defined in (35) and (52), 
and the on-axis dose distribution is given by (47) and (48). 

Figure 9 shows the relative lateral dose distribution at the Bragg peak (30.49 mm tissue 
depth) for a 60 MeV proton beam on muscle tissue. The initial beam width was taken to 
be 10 mm and the dose distribution was computed using (54). This result is in reasonable 
agreement was the results presented by NG, which show the lateral distribution in terms 
of the number of protons. Our results are broader in the wings of the lateral distribution. 
The contribution to absorbed dose in the lateral direction in Figure 9 approaches zero at 
approximately +1-9 mm. The distribution shown by NG approach zero in the lateral 
distribution at a lateral displacement of about +/- 7 mm. This result is to be expected 
since our radial broadening transport integral (i.e., (35)) tends to be larger than the TRIM 
results in the region between intermediate depths and the end-of-range, as discussed in 
section 4. 


Figure 10 shows the lateral distribution of relative isodose contours for the 60 MeV 
proton beam on muscle tissue. Broadening due to multiple scattering beyond the initial 
beam width is clearly evident at 26 mm tissue depth. As the end-of-range is approached, 
the beam rapidly broadens as the absorbed dose rapidly increases, and there is a well- 
defined cutoff in absorbed dose beyond the Bragg peak. 

7. Summary 

We have presented a 3-D formulation of ion beam multiple Coulomb scattering based on 
deterministic transport methods. The formulation includes an important step towards a 
self-consistent approach to coupling multiple scattering with ionization energy loss and 
straggling. Energy loss is coupled to lateral beam broadening due to multiple scattering 
by solving the Fokker-Planck transport equation with a depth-dependent, or equivalently, 
with an energy-dependent diffusion coefficient. The depth-dependent diffusion 
coefficient is evaluated by using the well-know range-energy relation. The range-energy 
relations itself is affected by multiple scattering processes as the actual particle path- 
length is greater then the penetration depth, which is accounted for in our formulation 
using the path-length correction to calculate an effective 1-D penetration depth. The 
methodologies and computation procedures presented in this paper will soon be 
integrated into the new 3-D version of HZETRN, which accommodates both space 
environment and ground-based laboratory boundary conditions. 
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We have tested our approach to ion beam broadening by calculating various broadening 
characteristics for a 60 MeV proton beam on muscle tissue. This scenario has important 
applications in radiation biology and ion beam therapy. Furthermore, Noshad and 
Givechi (2005) have calculated various broadening characteristics for the 60 MeV proton 
beam on muscle tissue using the TRIM Monte Carlo code from which we can compare. 
The results between our deterministic code and the TRIM simulations generally compare 
quite favorably for the 60 MeV proton beam on muscle tissue. Future efforts will be to 
extend our formulation to include inhomogeneous target materials and to incorporate 
beam attenuation due to nuclear absorption process, as well as a thorough transport code 
verification and validation effort. 
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60 MeV Proton Beam on Muscle Tissue 



Figure 1: Radial broadening of an initial infinitesimal 60 MeV proton beam on muscle tissue. The 
radial broadening is calculated out to 31 mm in tissue depth. The Bragg peak is located at 30.49 mm 
(see Figure 8). 
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Figure 2: Radial flux distribution for an initial infinitesimal 60 MeV proton beam on muscle tissue. 
The radial flux is shown along one of the lateral axis (i.e., either x- or y-axis). The radial flux is shown 
at various depths in the muscle tissue. 
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60 MeV Proton Beam on Muscle Tissue 
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Figure 3: Radial broadening of a 60 MeV proton beam on muscle tissue with an initial beam 
diameter of 2 mm. The radial broadening is calculated out to 31 mm in tissue depth. The Bragg peak 
is located at 30.49 mm (see Figure 8). 
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60 MeV Proton Beam on Musale Tissue 
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Figure 4: Radial flux distribution for a 60 MeV proton beam on muscle tissue. The initial beam 
diameter is 2 mm. The radial flux is shown along one of the lateral axis (i.e., either the x- or y-axis). 
The radial flux distribution is shown at various depths in the muscle tissue. 
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60 MeV Proton Beam on Muscle Tissue 
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Figure 5: Percent increase in the effective 1-D path-length due to multiple scattering for a 60 MeV 
proton beam on muscle tissue. Exiting ion beam assumed to have zero angular displacement, 
irrespective of lateral displacement. The percent increase in the effective 1-D path-length is 
calculated out to 31 mm. The Bragg peak is located at 30.49 mm (see Figure 8). 
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60 MeV Proton Beam on Muscle Tissue 
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Figure 6: Mean energy as a function of depth for a 60 MeV proton beam on muscle tissue. The mean 
energy is shown with and without the multiple scattering path-length correction applied. The Bragg 
peak is located at 30.49 mm when scattering effects are included. The Bragg peak is located at 30.73 
mm when scattering is not included. 
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60 MeV Proton Beam on Muscle Tissue 
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Figure 7: Energy straggling width as a function of depth for a 60 MeV proton beam on muscle tissue. 
Straggling is shown with and without the multiple scattering path-length correction applied. The 
straggling width is calculated out to 31 mm. The Bragg peak is located at 30.49 mm (see Figure 8). 
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60 MeV Proton Beam on Musale Tissue 



Figure 8: On-axis, relative dose-depth curve for a 60 MeV proton beam on muscle tissue. Path- 
length corrections due to multiple scattering are included (see (47)-(48) and the subsequent 
discussion). 
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60 MeV Proton Beam on Musale Tissue 



Figure 9: Lateral (relative) absorbed dose distribution at the Bragg peak (30.49 mm tissue depth) for 
a 60 MeV proton beam on muscle tissue. The initial beam width is 10 mm. 
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Relotive(%) laodose Curves: 60 MeV Preton Beam on Muscle Tissue 



Figure 10: Relative isodose curves for a 60 MeV proton beam on muscle tissue. The initial beam 
width is 10 mm. 
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